Modellization of hydraulic fracturing of porous 
materials^] 

F. Tzschichholz 1 & M. Wangen 2 

1 Department of Physics, Norwegian University of Science 
and Technology, N-7034 Trondheim, Norway 

EMail: frank@ical . uni-stuttgart. de 

2 Institute for Energy Technology, N-2007 Kjeller, Norway 
EMail: magnus@ife.no 

Abstract 

We review micro-structural fracture growth models suitable for the study of hydraulic 
fracture processes in disordered porous materials and present some basic results. It is 
shown that microstructural models exhibit certain similarities to corresponding theories 
of continua. These similarities are most easily demonstrated for simple crack geometries, 
i.e., straight cracks (finite size scalings). However, there exist even scaling relations which 
are completely independent of the particular employed crack structure. Furthermore it is 
demonstrated that disorder in cohesional/flow properties can influence the crack growth 
and the resulting fracture geometry in an essential way. 

1 Introduction 

The classification and prediction of breaking processes within solid materials 
forms an extremely difficult problem in engineering science as well as in solid state 
physics. This is even true for brittle materials exhibiting the simplest rheological 
behavior, i.e. no remanent deformations. Continuum mechanical treatments of 
brittle solids based on Griffith's theoryEJ are quite widespread and accepeted in 
scientific papers, and have been relatively successful for practical and engineer- 
ing purposes.0 There are, however, a couple of problems being associated with 
classical fracture mechanics limiting its range of application. Certain restric- 
tions arise from the assumption that the deformations can be conducted in a 
thermodynamical reversible manner. Connected to this is the basic assumption 
that the involved fracture processes are thermodynamically reversible too. This 

1 Reference: Advances in Fracture Mechanics, Fracture of Rock ed. M. H. Alibadi, (WIT 
Press, Southhampton, 1999), p. 227-61. 



is usually expressed by minimizing thermodynamic potentials with respect to 
crack/fracture growth, which is a perturbational scheme. The typical modelling 
situation is such that some governing field equations are describing the global 
state of the deformation (i.e. the free energy density) according to some bound- 
ary/loading conditions while the free surface energy field (intrinsic strength) is 
classically assumed to be a constant property throughout the solid material. The 
assumption of a constant intrinsic strength constitutes a quite strong idealiza- 
tion for most experimental situations because it eliminates the field nature of 
cohesional properties. The assumption of a constant theoretical strength is most 
easily justified for homogeneous materials, i.e., perfect single crystals. However, 
even for such classes of materials there exist certain limitations: a) the crys- 
tal's cohesion is in general non-isotropic (variation in crystal strength according 
crystallographic directions) and b) the cohesion is subject of thermal influences 
(lattice vibrations) introducing thermal fluctuations into the problem. 

The survey becomes more difficult for most solid materials of practical inter- 
est as they are frequently composite materials. Variations of cohesional properties 
occur already on macroscopic length scales comparable to those of the microstruc- 
tural elements. There is no simple argument to replace the theoretical strengths 
of the isolated solid phases by some average effective strength in connection with 
fracture problems. This is because the fracture and its associated natural ge- 
ometry advances on a microstructural level where the heterogeneity in grain ori- 
entation and strength is present. We mention this in order to illustrate that 
the assumption of constant cohesion/binding properties throughout a solid and 
therefore the representation employing a homogeneous continuum needs careful 
consideration. 

The investigation of geometrical properties of irreversibly grown cracks in the 
presence of disorder has recently received considerable attention amoung physi- 
cists. Basic ideas go back to statistical physics, theory of phase transitions, scaling 
theory, selfsimilarity and self-organized criticality. These concepts are described 
in good standard books of condensed matter physics, i.e. Chaikin & Lubensky.S 
Perhaps one of the most well studied phenomena exhibiting self-organized crit- 
icality is Diffusion Limited Aggregation (DLA), theoretically first implemented 
by Witten & Sander .□ Its outline is quite simple and it is well suited to exemplify 
geometrical selfsimilarity. 

Particles being small enough to be subject to Brownian motion diffuse in a large 
gas or fluid reservoir. The particle concentration is assumed to be very low (in- 
finite dilution), i.e., there are at each time step precisely one diffusing particle 
and one cluster where the diffusing particle is allowed to be attached to. If the 
diffusing particle touches the particle-cluster it sticks to it forever and a new 
diffusing particle is started far away from the cluster. The process is called diffu- 
sion limited because the particle diffusion is the rate limiting step in this reaction 
scheme. Figure [l] shows such a DLA cluster. The cluster has a dendritic structure 



Figure 1: Selfsimilar DLA-Cluster in two dimensions consisting of 10 4 particles. 
The cluster grows outwards 'absorbing' diffusing particles coming from remote 
distances. The fractal dimension is df = 1.71. With courtesy of S. Schwarzer 
(numerical simulation) . 

exhibiting branches on all length scales between particle and cluster length scale. 
This supposes that the average number N of particles belonging to the cluster 
follows a power law in cluster radius R, i.e., 

N~R d f. (1) 

The exponent df is called the fractal dimension of the cluster. The value for df is 
in general always smaller equal the embedding dimension; in the above example 
df = 1.71 < d = 2. The physical meaning of a power law relation Eq.(|l|) is 
that the structure does not contain any natural length scale and that it looks 
similar to itself (selfsimilar) under a change of scale. The reader interested into 
concepts and applications of Fractals is referred to References ||| and || Nowadays 
it is widely believed that the value of the fractal dimension does only depend on 
the symmetries of the governing equations but not on particularities such as the 
particle shape or the employed discretisation scheme in numerical simulations. 
Deviations from scaling behavior are anticipated for systems of finite size. 

The main interest in the DLA problem results from the fact that the gov- 
erning equation is the Laplace equation being one of the most widespread and 
simplest equations in physics. A diffusing particle sticks with probability propor- 
tional to dP/dn to the clustersurface with P being the solution of the diffusion 
equation. The sticking-probability is highest near the clusters extremities, so that 



Figure 2: Lichtenberg figure in an epoxy plate. The figure shows the electrical 
discharge pattern towards a notch at the right side of the plate.i 



the cluster is growing at a higher rate at the tips than in its interior. The tips 
screen the cluster interior to a certain amount leading to a diluted dendritic struc- 
ture. Other examples belonging to the same universality class as DLA are the 
dielectric breakdown (DBM) (see Fig. H\. the viscous fingering (VF) (see Fig. |a) 
and the Stefan problem of solidification. 3! All this problems lead experimentally 
and numerically to similar geometrical patterns. 

There has been considerable discussion whether fracture processes in presence 
of fluctuations would lead to DLA-like breakdown pattern too. The conducted 
investigations do, however, not support this idea in an unambiguous way. Rather 
the structure of the basic elastic equations (Lame or Cosserat equations) cannot 
be reduced to a Laplace type equation even in limiting cases for the elastic con- 
stants. The experiments of Van Damme and collaborators on twodimensionally 
confined clay pastes display a sharp transition between DLA morphologies (vis- 
cous fingering, compare Fig. ||a) and the much more spiky and diluted hydraulic 
fractures, e.g. see Fig. [|b. 

Further below we will present numerically obtained results mainly concerned 
with the problem of hydraulic fracturing of disordered materials in two dimen- 
sions. 

We will first consider the case in which the model material contains only elastic 
interactions. We will investigate the cases of a constant injection (borhole-) pres- 
sure and of a constant fluid influx. 

In a later chapter we consider a liquid saturated porous elastic material being 
fractured in a similar way. However, in such a case is the pressure distribution 
originating from the pore fluid flow directly coupled to the elastic equations. 



Figure 3: Fractal structures obtained after injecting water into Hele-Shaw cells 
containing clay pastes, (a) Viscous Fingering in a thin liquid paste (solid/liquid 
= 0.08). The fractal dimension is df = 1.65. (b) Hydraulic Fracturing in a thick 
paste (solid/liquid = 0.2). The fractal dimension is df = 1.43. Note that the 
fracture structure is more spiky than the viscous finger pattern. Both patterns 
were formed at the same constant injection pressurei!2l 

2 The Model 

In the following we will outline the employed model. We give a brief descrip- 
tion of the basic elastic and flow equations. Thereafter we explain in detail the 
employed boundary conditions. After this we demonstrate how heterogeneous 
cohesive properties are taken into account. Finally we present the breaking rules 
which contain the physics of the considered breaking process. 



2.1 Pore space 

There are at least two different ways of implementing porosity numerically. In 
both cases certain microstructural elements (beams, plates or shells) will form 
a stress carrying elastic backbone. One can think, however, of two different re- 
alizations for the microstructural elements. The elements could be 'made' of 
a material which already follows the constitutive laws for porous materials, e.g. 
Biot's theoryJHI This would be most satisfying from a continuum mechanical point 
of view as long as no fracturing occurs. However, if fracturing at the pore level is 
to be taken into account the above mentioned approach does not suffice (or even 
becomes inconsistent) and a micromechanical approach must be employed. The 
pore space/volume has to be modelled explicitely in that case. We have focussed 
on the latter implementation: the microstructural elements follow the classical 



constitutive equations of linear elasticity, and the pore space is viewed as the 
remaining vacant volume. The breaking of microstructural elements is associated 
with a coallescence of pore volume under question. Because the pore volume is 
modelled explicitely the couplings between elastic and flow interactions ought to 
be derived explicitely from the microstructural configuration. 

In this paper we consider a two-dimensional model-microstructure laying on 
a square lattice with properties of the pore space (hydrostatic pressure field) de- 
fined on its dual lattice.E3 



2.2 Elastic equations 

We consider the beam model (as defined in p. 232 of Ref. [i~3| ) on a two di- 
mensional square lattice of linear size L. This model is preferable compared 
to Hookean spring models from a microstructural point of view because elas- 
tic beams due to their finite crossection can be physically interpreted as fibers, 
while Hookean springs (having zero cross section) can not. In the two-dimensional 
beam model each of the lattice sites i carries three real variables: the two transla- 
tional displacements Xi and yt and a rotational angle ipi in the xy-plane. Neigh- 
bouring sites are rigidly connected by elastic beams of length I. The beams 
are assumed to have the same quadratic cross section, A = d 2 , and the same 
elastic behavior, governed by three material and geometry dependent constants 
a = I /(Ed 2 ), b = l/(Gd 2 ) and c = Ul 3 /(Ed 4 ) with E and G being the Young 
and shear moduli. With characteristic values for rock materials E = 10 11 Pa, 
G = 4 • 10 10 Pa, I = 10 -5 m and d = 10 -6 m the above three constants become 
a = l(T 4 (Pa m)~ l , b = 2.5 • 10~ 4 (Pa m) _1 and c = 0.12 (Pa m)" 1 . When a 
site is rotated (<fi ^ 0) the beams bent accordingly always forming tangentially 
90° angles with each other. In this way local momenta are taken into account. 
For a horizontal beam between sites i and j one has the longitudinal force acting 



at site j, 




(2) 



the shear force 




(3) 



and the flexural torque at site j 





using a = 1/a, f3 = 1/(6 + c/12) and 5 = f3(b/c + 1/3). The corresponding equa- 
tions for vertical beams are similar. The outline of the beam network is displayed 
in Fig.||. We will not consider inertial or gravity forces, however, the flow of a 
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Figure 4: Schematic representation of the beam model (left) rotation of one site, 
(right) beam flexed due to the angles at its extremities. 

pore fluid amounts in general to non-vanishing body forces. If there is a pressure 
drop p2 — pi across a beam (pore wall) the resulting line force density can be 
replaced by a shear force dipole, S = (p± — p2)ld/2, acting at both ends of that 
beam. This follows directly from the equilibrium conditions for beams. Hence 
the dipoles couple the pressure field in form of body forces to the mechanical 
equilibrium conditions. 

In mechanical equilibrium the sum over all elastic and body forces (torques) 
acting on site j must vanish. Thus if the pressure field is known one can calculate 
the displacements. 

To determine the hydrostatic pressure field one needs to consider flow equa- 
tions. 



2.3 Flow equations 

The beams exhibit a typical length I and thickness d and the liquid saturated 
equilibrium pore volume is simply [l — d)d 2 . This states explicitly the considered 
material to be a binary system, i.e. any volume element is either occupied by 
elastic material or by pore fluid. By varying the ratio of I jd different values for the 
volumetric porosity can be achieved (ranging from nearly zero to nearly one). In 
order to obtain most simplified flow equations, we make the following assumtions: 

a) the pore fluid is incompressible (which appears well justified for water) and 

b) the fluid velocity relative to the bulk is proportional to the gradient of the 
pressure field (Darcy's relation). The continuity equations for solid and liquid 



mass density lead to the required pressure equation (in the continuum limit); 



div^— grad pj = — div u, (5) 

with p being the pressure, u the displacement vector, k the mechanical perme- 
ability and \x the fluid viscosity. This equation is a Poisson equation for the 
pressure field for constant k and /i. The right hand side of Eq.(|5|) represents in 
general a time-dependent source term for the pressure field. It characterizes the 
local rate of relative volume change due to elastic deformations, and hence the 
net fluid flux itself. In general these net fluid fluxes will control the pressure 
field's time-evolution. However, in this work we will only consider the stationary 
case. The physical motivation for this is the observation, that the characteristic 
internal relaxation time r = [i/E for water-rock systems is typically of order 
10" 6 s. It appears justified to assume that the time periods between sucessive 
breaking events are much larger than the internal relaxation time r, and hence 
that breaking events are triggered by stationary pressure and displacement dis- 
tributions. In connection with Eq.([|) this implies to balance the net fluid flux 
for a given cell i to zero (Laplacian); 

— kij(pi — Pj) + boundary conditions = 0. (6) 
^ j 

Here ku denotes the local permeabilities across the beams (pore walls). 

The pressure equation Eq.(|6|) is solved without referring to any mechanical 
properties, which is directly a consequence of explicitely dropping all dynamic 
terms. Its solution for given boundary conditions represents the body force under 
which the elastic equations, among their own boundary conditions, are solved. 
Formally the correspondingcontinuum problem is analogous ±o the standard 
problem of thermoelasticitylia which was already noted by BiotEJ. From the re- 
sulting stresses/strains and breaking thresholds (see Sec.|2.5|) we then determine 
which element(s) are broken next. 

Here we would just like to mention what happens to the pressure equations, 
when a beam is broken. Because the beams represent the pore walls, and because 
the fluid is assumed to be incompressible the pressure within a connected crack 
has to be constant. In this sense a crack is a macro-pore. Numerically this can 
be accomplished in different ways. We have chosen to set the permeabilities for 
broken beams very high, i.e. kij = 10 9 . Hence the pressure drop along a crack is 
negligible small. 



2.4 Boundary conditions 



The results pre sented in this review are mainly for hydraulic fracturing conditions 
in impermeablell3ll^l and permeableS materials. The only exception is Sec. || 
where we compare results for tensile and hydraulic loading configurations .t3 The 
corresponding tensile boundary conditions are therefore discussed in Sec. || 

Let us make few remarks about terminology. A boundary will be called 
'free' if its normal stress component vanishes and will be called (hydrostatically) 
'loaded' if the normal stress component follows some known non-zero stress dis- 
tribution (von Neuman condition). A boundary is called 'fixed' or 'clamped' 
if the displacement values are assigned to the boundary (Dirichlet condition). 
A particularity in calculations are the sometimes employed 'periodic' boundary 
conditions for displacements /stresses. They enforce periodic displacement /stress 
field solutions. Periodic boundary conditions are especially useful in numerical 
calculations when the asymptotic behavior for infinite systems is of interest. We 
will encounter them in Sec. ||. 

For most physical fracture problems the boundary conditions for the displace- 
ments are von Neuman conditions, or mixed boundary conditions. In presence of 
an additional flow equation boundary conditions for the fluid pressure field need 
to be specified. In our case homogeneous Dirichlet conditions are sufficient, i.e., 
the hydraulic crack pressure, p, and an equilibrium pressure, po, for the external 
boundary. 

In Sec. 3^2 and Sec. || we describe simulations where the hydraulic crack is 



'loaded' and propagated by a constant crack pressure. The crack pressure does 
neither vary in space nor in time. The considered external boundaries (lattice 
edges) are either free or periodic. 

A, for practical purposes, somewhat more realistic boundary condition prop- 
agating a hydraulic crack is a constant fluid injection rate, qi n . To keep things 
tractable we only consider an incompressible hydraulic fluid. The boundary con- 
dition for an impermeable crack, Sec. [O], is a particular case of the more general 
permeable crack boundary condition, see below. A constant fluid injection rate 
qi n corresponds to a certain driving pressure difference, p — po, where p denotes 
the hydraulic pressure within the crack and po corresponds to an equilibrium 
pressure. We have used po = 0. The hydraulic pressure does depend on qi n and 
the crack opening volume. For permeable cracks/materials p also depends on the 
hydraulic losses. In the following we will outline how to determine the hydraulic 
crack pressure p for given injection rate qi n and crack structure (contour) d V. 

Because for an incompressible fluid the conservation of fluid mass directly 
implies a conservation of fluid volume, it is possible to formulate a continuity 
equation for the crack opening volume V, 

dV f 

— + / v-dA = q in . (7 
at Jdv 



Here v denotes the Darcy field integrated over the crack surface d V with surface 
element dA. As source term appears the rate of injected (fluid) volume, q% n . 
Because the Darcy vectors are linear funtionals of the pressure, v = —£;///• gradp, 
and the crack opening volume is a linear fuctional of displacements, V = Jq V u ■ 
dA, and because all considered equations are linear, a scaling of the form p = Xp* 
implies v = Av*, u = Au* and V = XV*. Note that such a scaling is in general 
impossible if one considers Eq. (||) instead of the Laplacian, Ap = 0, as flow 
equation. 

There exists an interesting application of Eq.(^) for undercritical, non-growing 
cracks. With V* and v* depending only on the constant crack contour d V and 
crack pressure p* one obtains a first order differential equation for the scaling 
factor X(t), 

dX 
~di 

with general solution, 



, , u*-dA + X- v*-dA = q in , (8) 
at JdV JdV 



X(t) = Aoo + (A - Aoo) • er*/ T , (9) 



using 



> Qin . _Po J dv u*-dA 

A °°- / av V-dA' A °"^' T ~ J dv v*.dA- (10) 

The important point is that one can choose an arbitrary boundary value p* 
within a calculation; the resulting p = X(t) ■ p* will be unique. However, one 
needs to calculate the contour integrals in Eq. ( JlOj ) in order to do the rescaling. 

The magnitude Aoo does depend only on the injected fluid flux qi n and the 
total fluid flux through the crack contour. It does not depend on elastic properties. 
The typical time scale r, above which the crack pressure and the crack opening 
volume tend to approach their asymptotic values, is given by the ratio of crack 
volume and fluid volume flux through the crack surface. It does not depend on 
the injected fluid flux. Again we would like to stress that the crack contour is 
assumed to be constant. However, Eq.(||) is quite general in the sense that all 
dependencies on the actual crack contour have been implicitly taken into account 
within the constants Aoo an d r. Equation @ should therefore hold for arbitrarily 
shaped cracks. 

There are two limiting cases for the practical important situation Ao ~ 0. 
For large times, 4 > t, the scaling factor X(t) — > A^o becomes constant and so 
do the crack pressure Poo — QinT ' (p*/V^*) and the crack volume Voo — QinT- 

The 

other limiting case is for small times, t <§C r, where the scaling factor becomes 
linear in time, X(t) = ^ • t with corresponding pressure p = qi n y^ ■ t and 
volume V = qi n ■ t. The latter case also holds for impermeable materials because 
lim fc ^o = qin/V* ■ t. 



Equations (0) and @ hold for arbitrary crack contours, i.e., branched or 
ramified cracks. 

Equation fj] has a corresponding discrete update scheme. One just needs to 
be careful with the time derivative in Eq. ([?]) because the scaling factors A are in 
general different before and after the time increment, At, 

\(t + A t) = X(t) ■ (l - ^A t) + f^A t. (11) 



V* > V 

Here D* = JgyV* ■ dA and V* = f dv u* ■ dA denote the crack contour inte- 
grals/sums for the Darcy flows and displacements respectively, as obtained from 
the solution of the boundary value problem at time t with boundary value p* 7^ 0. 
As soon as the crack growths and its boundary changes one needs to recaculate 
the flow and the elastic equations in order to obtain the corresponding D* and 
V*. The ratio D* /V* being the inverse of a typical time will therefore change 
during the crack growth process. 
We summarize: 

1. Choose an initial crack boundary. 

2. Set the initial value A(0), i.e., A(0) = (closed crack). 

3. Solve the pressure equation Eq. @ for an arbitrary hydraulic pressure p* 
(boundary value for the current crack boundary). 

4. Insert the corresponding pressure gradients into the elastic equations Eqs. (§)- 

as inhomogeneities and solve the elastic equations. 

5. Calculate from the obtained elastic solution the stress/force distribution 
a*. 

6. Calculate from the elastic solution the crack opening volume V* and from 
the pressure solution the total Darcy flux D*. 



7. Update A according to Eq.(|ll|). 



Use the updated stresses/forces, a = Xa*, in the breaking rule (see Sec. 2.6 ) 
and break the corresponding element(s). 



9. Repeat steps pfjq until a stopping criterion is fulfilled. 



2.5 Heterogeneities 

The problem of fracture processes in elastic media has undergone much attention 
in the last years, in particular the role of disorder in such processes. @H@ An 
important question of interest is how a variation in the local properties affects 
the global elastic properties in comparison to ideal homogeneous media. It was 



recognized that elastic properties could considerably be changed and therefore 
interest was attracted how failure mechanisms like fracture, rupture or damage 
processes are influenced, changed or possibly controlled by disorder. 

In the literature are two essentially different ways to take the disorder into 
account .Ej 

The first possibility is to use exclusively the strain / stress field in order to cal- 
culate the probability for breaking a local constitutive element. This is known as 
stochastic disorder. Prominent scalar examples for this kind of disorder are known 
in the literature as Diffusion Limited Aggregation (DLA), Dielectric Breakdown 
or Viscous Fingering.QO We will present results for stochastic disorder in Sec. ||. 

A second different kind of disorder is to initially assign a statically set of 
random values to certain material specific constants (structural disorder). Once 
these constants like strength, bending thresholds, elastic properties and others 
have been chosen the system behaves in a fully deterministic way. Corresponding 
results are presented in Sec. |5[ 



2.6 Breaking rules 

Lattice fracture models like the beam model are spatially discrete. This is be- 
cause one would like to have direct access to each constitutive elastic element. 
The medium is represented by a lattice and nearest neighbouring lattice points 
are connected by constitutive elastic elements called bonds. The bonds represent 
interactions not on an atomistic scale but rather on a mesoscopic scale, for exam- 
ple the elastic coupling between two adjacent grains where the grains themselves 
are linear elastic bodies. New features are the definitions of physically motivated 
"selection" and "breaking" rules for lattice bonds. Either a bond is broken ac- 
cording to some breaking rule and stays broken for the future or it remains intact 
and is further susceptible for the breaking process. The simulation ends when the 
elastic lattice fully breaks apart. Moreover, one has to define in the selection rule 
those bonds which could break in principle at a given stage of rupture. In single 
crack models only the surface bonds of the crack are allowed to break. However, 
in the many crack version in general all bonds are eligible. This allows in principle 
the nucleation of several microcracks in the system. We have used for our simu- 
lations the first kind of selection rule. In a recent published workE-3 it was shown 
that if the disorder for the breaking thresholds is not to weak there are two dif- 
ferent regimes in the macroscopic breaking characteristic. For small macroscopic 
strains first the weakest bonds are broken without leading to catastrophic failure. 
The breaking process is stable in this region and one has to increase the external 
applied stress in order to proceed breaking bonds. In finite systems however, this 
disorder controlled regime must end at the breaking point. The breaking process 
localizes, becomes unstable and a single large crack controlls the entire breaking. 
It was shown that for different probability distributions of breaking thresholds 



and different loading conditions two nontrivial possibly universal exponents are 
sufficient to rescale the macroscopic breakdown characteristics with respect to 
the system size. Breaking rules which determine the breaking conditions are cru- 
cial in lattice models because one has to put additional information about the 
fracture mechanism into the model which could not be obtained from elasticity 
theory. The employed breaking rules are described in the corresponding sections 
Sec. ||, Sec. and Sec. ||. 



3 Fracturing without disorder 

We investigate both the case of uniaxial loading due to remote tensile forces and 
the loading due to a hydrostatic pressure acting from the interior on the crack 
surface. The former situation is typically encountered when a material is pulled 
at its outer faces, whilethe latter one can be found in engineering applications 
of hydraulic fracturing. c3 We will neglect inertial forces, assuming a sufficiently 
slow fracture process. The cohesive properties of the beams are considered to be 
homogeneous, i.e. the beams behave linear-elastic up to their theoretical strength 
F t h (cohesion force), above which they irreversibly break (zero elastic constants). 
At the beginning of each simulation we break one vertical beam located at the 
center of the lattice. This is the initial crack. Corresponding to the employed 
boundary conditions described below we calculate the internal force distribution 
Fij of beams connecting sites i and j using a conjugate gradient method.E3 From 
this longitudinal forces we determine and break the most over stressed beam, i.e. 
the beam carrying the highest force, Fi j = max/j^i F^. The maximum force 
allows us to calculate a scaling factor A = F t h/Fi o j = ath/&i jo which in turn 
is used to determine the macroscopic breaking stress a c or breaking pressure P c 
(measured in units of oth) necessary to fulfill the local fracture criterion (Ti j = 
(Tth- Such scaling is possible due to the linearity of the employed elastic equations. 
Breaking the beam connecting sites iq and jo destroys the balance of forces at 
those sites and one has to relax the system to its new equilibrium configuration. 
Then the above steps are repeated until the lattice breaks apart. 

It is well known that the highest stress enhancement factors occur at the 
crack tips for the pressure as well as for the tensile problem. Because we consider 
only homogeneous cohesion forces the terms 'most over stressed' and 'highest 
stress' are equivalent and the cracks grow linear in direction of the highest local 
tensile force. In our case the direction of crack growth is the x-direction because 
of the uniaxiality of imposed boundary conditions, see below. 



3.1 Tensile fracturing 

To impose an external strain we attach at the bottom of the lattice a zeroth 
line on which for all sites Xj = yi = <pi = are fixed, and on the top we attach 



a (L + l)st line on which all sites have the same fixed values Xi = <pi = and 
Hi = 1. With this boundary conditions the external displacement in y-direction is 
fixed to unity (Dirichlet boundary conditions) and one can imagine them as being 
represented by rigid bars attached at the top and the bottom of the lattice. For 
the boundary conditions in x-direction we consider periodic as well as free lattice 
boundaries. The externally applied stress, do, necessary to maintain a unity 
displacement between the two rigid bars can be easily calculated by summing up 
all yi displacements over the first line, <7o = -jfa J2 Vii where a is the longitudinal 
force constant, A the cross section of the beam and L the (dimensionless) lattice 
size. Employing the above defined scaling factor A we express the externally 
applied breaking stress, a c , in the form, 

o c = Acr , (12) 

or in dimensionless form, 

— = -J-T 1 Evi. (13) 
a th F iojo L 

With this notation the external breaking stress a c is just the stress necessary to 
yield the local breaking condition (Ji j = ath- 

We will use Eq.(^) for our finite-size scaling purposes. It expresses the 
breaking stress, measured in units of cohesion stress, in terms of the beam force 
Fi j at the 'crack tip' and the average force ^ J2 yt on a beam due to the global 
unit displacement. In general both terms depend on the number of broken beams 
N (crack size) and on the lattice size L, see Sec. |3.3.1 . 

3.2 Hydraulic fracturing 

Experimentally hydraulic fracturing is encountered when an incompressible fluid 
is injected under high pressure into an existing crack. The characteristics of 
this particular loading mode is that the loading of the crack only happens on the 
crack surface itself, but not on remote surfaces like in the case of tensile fracturing. 
Recently hydraulic fracturing has been numerically addressed by introducing for 
each broken beam a force dipole of strength Fq simulating an inner pressure Po = 
.Fq/AEJ'ES In the present section we will follow this implementation. The basic 
assumptions are a) the pore fluid is highly compressible (gas), b) the permeability 
is negligible and c) the injected fluid is incompressible. Starting from one vertical 
broken beam (zero elastic constants, one vertical force dipole) we calculate the 
internal stress distribution Fij using the conjugate gradient. Calculating the 
scaling factor A as described above we determine and break the beam carrying 
the highest stress, Fi Q j Q . Thereafter the lattice has to be relaxed again (now two 
broken beams, two force dipoles) and one repeats the above mentioned steps until 
the lattice breaks apart. The crack growth is single connected. As always the 



two vertical beams in front of both crack tips are most overstressed and broken. 
Hence the existence of a connected path for the liquid is garantued. 

We have performed calculations for free external boundaries in x- and y- 
direction as well as for periodic boundaries in x-direction and free boundaries in 
y-direction. As we only impose stresses the elastic solution (displacement field) is 
unique up to a general translation and rotation. Without restriction we therefore 



fix the displacement of an arbitrary lattice site to zero. Analogous to Eqs. (12) 
and ( |l3| ) one has for the breaking pressure, 

P c = AP , (14) 

and 

respectively. It should be noted that <Jth is the tensile strength and not the 
compressional strength. In fact the stresses at the tips are tensile stresses as the 
crack is hydraulically opened. 

We will use Eq. ( |T5| ) for the finite-size scaling of the hydraulic problem. 

3.3 Results 



As explained in Sec. 3.1 and 3.2 we obtain from our simulations the macroscopic 
breaking stresses and pressures, a c (N, L) and P C (N, L), measured in units of the 
theoretical strength ath, in terms of the number of broken beams N and the 
employed lattice sizes L. As we are interested in the finite-size scaling properties 
of the breaking stresses (E c = a c ) and pressures (E c = P c ) we introduce two new 
functions, *(iV) and $(N/L), 

X c {N,L) = a th V{N)<f>(N/L). (16) 

The function $(N/L) describes the finite-size corrections to the scaling \P(iV) 
between the breaking stress/pressure and the crack length for given boundary 
conditions in the infinite system. The quantity N/L is for our purposes the most 
suitable scaling variable as it represents the length of a linear crack measured in 
units of the employed lattice size. In general both functions, \£ and $ will de- 
pend on the employed boundary conditions. The dependence of the asymptotic 
scaling ^ on employed boundary conditions appears to be contrary to what is 
expected from the theory of critical phenomena. In fact the asymptotic scalings 
^ for breaking stresses and pressures of continua without internal length scale 
are identical. However, as we will see below this is no longer necessarily true for 
systems with internal length scale. One obtains rather finite-size corrections orig- 
inating from this scale which, depending on the employed boundary conditions, 
can result in very broad, non universal transient regions. 



In order to show numerically that the finite-size scaling function <&(N/L) 
exists, one needs to know (or to conjecture) the explicit form of *$>(N), i.e. the 
asymptotic behavior of S c as L — > oo. 

We will access the latter information from simple continuum mechanical con- 
siderations of symmetric elasticity. 



3.3.1 Tensile fracturing 

Previous continuum mechanical scaling considerations have successfully been ap- 
plied to fracture models. It has been shownlla that the breaking stress a c for 
tensile fracturing of a finite, periodic beamjattice follows very well the known 
'Tangent-Formula' of fracture mechanics c3'E3, 

_^-./y_ cot __. (17) 

By comparison with Eq. (|l6|) we see that the asymptotic scaling for the breaking 
stress o~ c is given by ~ N~ 1 / 2 , a result that originally was found by Griffith. 
In Fig. | we show the finite-size scalings for the numerically obtained tensile 
breaking stresses (a) for free boundary conditions and (b) for periodic boundary 
conditions. For both curves the data collapse is quite acceptable. While for 
the periodic problem the finite-size correction &(N/L) follows Eq. (|l7])EI, no 
analytical expression for is known in the case of free boundary conditions. 
However, the differences for the breaking stresses at given values of N/L are less 
than O.lut/j, see Fig. [5L and Eq. (|17|) represents a reasonably good approximation 
for finite solids for many practical purposes. 



3.3.2 Hydraulic fracturing 

In general analytical solutions of crack problems are difficult to obtain and even 
for two-dimensional problems in infinite, semi-infinite or periodic domains one is 
often faced with complicated integral equations. The reader interested in this is 
referred to MuskhelishviliEJ and SneddonEE 

It has been shown analytically that the breaking pressure P c of a semi- 
infinite two-dimensional continuum with periodic boundary conditions is given 
by Eq. (ff^) where one has to replace a c by -P c £3 We show in Fig. |6| the analogous 
plot to Fig. H for the breaking pressure. One does not obtain a data collapse, 
indicating that \& ~ iV -1 / 2 , see Eqs. ( |l6|) and (17), is not the scaling of the 



breaking pressure in the asymptotic limit of infinite lattices. 

This result is somewhat unexpected, because the continuum limit predicts 
identical scaling relations for breaking stresses and pressures and on the other 
hand we find good agreement between lattice and continuum scalings for the 
tensile problem. 
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Figure 5: Finite-size scaling for the rescaled breaking stress o c j 'athN 1 ^ 2 as a 
function of the rescaled crack length N/L, (a) for free boundaries in x-direction 
and (b) for periodic boundaries in x-direction. Lattice sizes: (A) for L = 20, (+) 
for L = 40 and (o) for L = lOoH 




However, the proposed model for hydraulic fracturing is defined on a lattice 
and not in a continuum. As lattices always show additional structure, higher 
order gradient terms of the continuum displacement field appear in the elastic 
solution, which have no counterparts in simple continua.0@E3 

In the following we present an argument why the lattice finite-size scalings 
for the pressure and tensile problem are different. 

In the tensile problem the loading forces are acting at the lattice boundaries 
remote to the stress free crack surfaces. Therefore local stresses close to the crack 
tips are carried through nearly the whole elastic volume. Contrary to this in the 
hydraulic problem loading only happens due to force dipoles acting at the crack 
surfaces close to the crack tips. One therefore might expect that the breaking 
pressures P c show a more pronounced '(micro structure) lattice behavior' than 
the breaking stresses a c - For very large cracks these differences should diminish 
as the discrete loading approaches a force density. 

In order to extent a linear crack of length N in an infinite continuum one has 
to impose a breaking pressure P c ~ N^ 1 / 2 . However, employing a single double 
force of magnitude -F^ne finds, F c ~ TV 1 / 2 for the critical loading force (related 
Boussinesq problem) .E3 

We argue now that for a single crack in an infinite disk loaded by equidistant 
force dipoles the asymptotic scaling for \E r is equivalent to that of an infinite 
lattice. Knowing ^ we plot the finite-size scalings for the investigated lattices. 




Figure 6: Plot of the rescaled breaking pressure P c / (JthN 1 / 2 as a function of 
the rescaled crack length N/L for periodic boundaries in x-direction. In contrast 
to Fig. H we do not find any data collapse, indicating that ^/(N) ~ iV -1 / 2 in 
Eq. (16) is not the appropriate asymptotic expression for the breaking pressure. 
Lattice sizes: (A) for L = 20, (+) for L = 40 and (o) for L = 100 H 



The asymptotic continuum scaling for the breaking pressure (as calculated 
from equidistant force dipoles) is most easily obtained using the complex stress 
function of Westergaard.E3'E3 We find for the asymptotic scaling in an infinite 



disk, 



with 



7T N 1 / 2 

•™ a immry (18> 

(2V-l)/2 



The sum f(N) converges rapidly towards arcsin( j^~j ) and for large ./V towards 
7r/2. Hence in the limit of a large number of force dipoles we obtain the afore- 
mentioned asymptotic Griffith scaling for the breaking pressure, P c /&th = -/V -1 / 2 . 
In Fig. |?] we show the resulting finite-size scalings for the breaking pressure, 



based on Eqs. ( jig) and (18). In this plot a reasonable scaling is obtained, which 
proves our earlier statement that the asymptotic limits (L — > oo) of the breaking 
pressures for continua and lattices are different. 

Comparing Fig. ^ with Fig. [| one also notes a significant deviation from the 
scaling form Eq. (|j~7D, demonstrating a modified correction to scaling behavior. 
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Figure 7: Finite-size scaling for the rescaled breaking pressure P c /oth^ \N) 



with ^ 1 {N) = 1 1+ ^ 1 /2 jV ^ , see Eq.(18), as a function of the rescaled crack length 
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4 Hydraulic Fracturing with stochastic disorder 

We briefly discussed in the introduction, Sec. |l], Diffusion Limited Aggregation 
and its associated fractal aggregates. Essentially DLA clusters growth by adding 
one particle per "simulation step" to the clusters surface with a probability, pij, 
proportional to the local field gradient |VP| obtained from Laplace's equation, 
A P = 0. Similar considerations apply for the quasistatic growth of electrical dis- 
charge patterns, where an elementary discharge event is assumed to happen with 
a probability proportional to the local electrical field.u Typically the boundary 
conditions are Dirichlet conditions for the external and the cluster surface. 

There have been attempts in literature to generalize the DLA model in sev- 
eral ways. A main proposition was the introduction of a growth exponent, 77, 
describing how strong the growth probabilities are correlated to the gradients, 
Pij ~ iVPp. For 7] — > 00 the probability distribution becomes arbitrarily narrow 
(a delta peak) and the model is completely deterministic, i.e., the cluster grows 
always at the surface location of highest field gradient (current). The resulting 
structures are usually linear in two dimensions. The case 77 = 1 corresponds to 
fractal DLA cluster growth, as already mentioned above. For rj = the clus- 
ter growth is completely controlled by fluctuations and independent of the field 




Figure 8: (Left figure): A typical crack obtained with a beam model for hydraulic 
fracturing at constant pressure, on a square lattice of 250 x 250 with periodic 
boundary conditions. All beams under tension break with the same probability. 
The crack consists of 2200 broken beams. The greyscale represents the difference 
of the diagonal elements of the hydrostatic stress field. (Right figure): Only the 
crack structure.il 



solution. The corresponding structures are usually spatially compact. 

It is important to note that the growth rule and in particular the employed 
growth exponent r\ is completely independent of the investigated field equation. 
Rather the growth rule bases on plausible considerations. 

There has been subsequently work been done to apply stochastic growth 
rules onto elastic systems, especially applying them to crack growth problems. 
The "Gradient" in the growth rule is then typically replaced by the hydrostatic 
strain (eigenvalue) component. For certain classes of tensile problems fractal 
crack growth was found on two dimens iqna ] Hookean spring networks, exhibiting 
a fractal dimension dj ~ 1.6 for rj = l.Or 2 !! 23 ! This value for df is relatively close 
to the one observed in a DLA problem. However, the first author of this review 
conducted simulations on two dimensional beam lattices (see Sec. ^2) where a 
small preexisting crack was loaded at constant pressure from inside, and found 
non- fractal rather compact crack growth, already for 77 = 1.0 It was concluded 
that compressed microstructural elements (beams) should not be considered for 
the breaking process because of the very asymmetry of interatomic potentials. 
With this modified breaking rule also non- fractal cracks were found for r/ = 1. 
The cracks were, however, now linear in shape. In order to obtain in a hydraulic 
fracture simulation persisting branching cracks, much stronger fluctuations need 
to occur than obtainable with 77 = 1. This is not so surprising as almost the 




Figure 9: (Left figure): Double- logarithmic plot of the number of broken beams 
N as a function of the average crack radius R for free boundary conditions. Note 
the power-law behavior, N ~ R f , indicating fractal crack growth. (Right figure): 
Plot of the successive slopes df of Ihe data, indicating an average value for the 
fractal dimension df = 1.56 ± 0.05.&3 



whole elastic bulk is in a compressed state, while just around the crack tips there 
exist very narrow cones of elements being in a stretched state. In Fig. || we show 
the crack of a single, very large simulation for the case of maximum fluctuations, 
i.e., r\ = using the above mentioned asymmetric breaking criterion. Note that 
the breaking probability does not depend on the strain magnitude, however, all 
compressed beams are assigned a breaking probability identically to zero. The 
hydraulic pressure inside the crack was kept constant and the external boundaries 
are periodic in displacements, i.e., spanned on a torus. Furthermore the system 
is assumed as impermeable and there are no losses of the hydraulic fluid. The 
resulting crack structures are scaleinvariant (fractal) as can be seen in Fig.^. The 
measured fractal dimension df = 1.56 ±0.05 is slightly below the numerical value 
for tensile problems, and slightly above the experimental value df = 1.41 obtained 
from hydraulic fracturing of clay pastes in two dimensions, compare Fig.||. 

The total crack opening volume of the ramified cracks has been calculated 
and monitored during the growth, and plotted versa the average crack radius. 



This is shown in Fig. 10. The crack opening volume for fractal cracks shows 
the same dependency on (average) cracks size in comparison to straight cracks, 
namely V ~ R 2 . 

We note that the foregoing considerations and results are for impermeable 
elastic materials, i.e., the influence of possible effects due to hydraulic fluid flow 
(hydraulic losses) have been neglected. 
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Figure 10: Plot of the crack opening volume V as a function of the average crack 
radius R for fractal cracks with dimension df = 1.56 using (A) free boundary 
conditions and (□) periodic boundary conditions. For comparison we show also 
results for a straight crack (o) employing free boundary conditions. In all three 
cases we find the behavior V ~ R 2 (in d = 2). El 



5 Hydraulic Fracturing with structural disorder 

The concept of a local cohesion strength has been used in a number of papers.0 
One assumes that a deformed elastic element connecting sites i and j breaks 
above a certain material specific threshold force F+l ('cohesion strength'). This 
kind of disorder is also called deterministic or quenched disorder in the literature, 
because once the breaking thresholds have been chosen the model behaves in a 
completely deterministic way, contrary to the model described in the foregoing 
section. However, we prefer the term 'structural disorder' as a spatial variation 
in material properties is frequently associated with certain micro-structural par- 
ticularities within a material. If the inner stresses are above their thresholds 
the beams are broken and removed, i.e. their elastic moduli are set to zero. 
Since the cohesional strength for compression is much higher than for tension 
we assume that compressed beams cannot break. El The breaking thresholds are 
initially distributed randomly according to some probability density function, 
p(F t f l )- At this point it would be ideal to employ empirical strength distributions 
of the considered material for p(Fth)- However, at the present stage of research 
we are basically interested in generic features of such crack growth models. The 



Figure 11: Typical hydraulic crack for strong disorder, r = —0.7 and (F t h) = 
0.01, on a lattice of size L = 150. The crack consists of 629 broken beams after 
1500 time steps using a hydraulic flux q^ n = 0.05. The injection point is the 
center of the lattice, We have used periodic boundary conditions in vertical and 
horizontal directions .c3 

width of the strength distribution is most easily controlled employing power-law 
distributions, 

p{F th ) ~ F[ h , (20) 

with F t h £ [0, .Fmax] an d r > —1. Negative exponents r are used to describe 
strong cohesive disorder (broad distributions) while large positive exponents cor- 
respond to weak disorder (narrow distributions). It is convenient to express the 
normalization factor and F max by the distribution's expectation value {Fth) an d 
the exponent r. This allows to investigate disorder effects due to deviations 
around {F t h) for constant {F t h)- 

In the following we will discuss the modelling of crack growth in imperme- 
able and permeable materials separately because the employed breaking criteria 
slightly differ for both cases. 

5.1 The impermeable medium 

The hydraulic fracturing of an impermeable solid represents a limiting case due to 
the vanishing hydraulic losses. For a porous material saturated by an incompress- 
ible fluid this would imply that only pure shear deformations would be possible. 
This follows directly from the time-dependent flow equation, i.e., for k — > or 
[i — > oo Eq. (H|) reduces to the condition divu = const, preserving the volume of 
a deformation. This is not a very typical experimental situation. However, if the 
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Figure 12: Log-log plot of the number N of broken beams versus the average 
radius R(N) for two different distributions of breaking thresholds: (o) r = —0.7, 
averages over 60 cracks, fractal dimension df = 1.44±0.10; (+) r = —0.5, averages 
over 53 cracks, fractal dimension df = 1.39 rfc 0.10. For all simulations we have 
used an average cohesion value (Fth) = 0.01 and a linear lattice size L = 150.&3 



pores are filled by very compressible gases one might even neglect the deformation 
dependent volume forces for impermeable solids. 

In the following we consider the case where an incompressible fluid hydrauli- 
cally propagates a crack through an impermeable material. We require that the 
pores are filled by a very compressible substance compared to the solids compress- 
ibility. Such situation is approximately encountered for magma driven cracks£3 
We consider a constant fluid injection rate, qi n = 5T0~ 2 Z 2 d, into the hydraulic 
crack and a vanishing pore pressure, po = 0. The boundary conditions have 
been discussed in Sec. |2.4j , i.e., see Eq. ( |il~l) for vanishing Darcy losses, D* = 
The simulation procedure follows the mentioned steps [l||9] in Sec. 2.4 with 



0. 

the exception that in the present case no pressure equations need to be solved. 
Instead the 'test pressure' p* does appear directly as inhomogeneity for the elastic 
equations along the crack contour. Initially one small crack (one broken vertical 
beam) is prepared on the grid center. For the distribution of breaking thresholds, 
Eq. (f20|), the values (Fth) = 10~ 2 (arbitrary units), r = —0.5 and r = —0.7 
(strong disorder), and r = +oo (no disorder) are used. After each time increment, 
At = 1, all beams carrying tensile forces larger than the associated breaking 
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Figure 13: Linear plot of the pressure P inside the crack versus time t. The 
record corresponds to the simulation displayed in Fig. 11. ll 



thresholds are broken, i.e., their elastic moduli are set to zero. 

For r = +oo (no disorder) straight cracks are observed andJjie hydraulic 
crack pressure drops in form of a power-law in time, P(t) ~ t^ 1 / 3 .^ 

For strong disorder the situation becomes more complicated. In Fig. [ll| we 
show a calculated sample crack. The obtained crack is a branching structure 
exhibiting many small scale branches but only a very few branches on the larger 
length scale. A comparison with the crack structure shown in Fig. || reveals a 



much lower crack density in the present case. This is directly measured in Fig. 12 



by_plotting the number of broken beams N in dependence of the average crack size 
i?,E3 yielding the cracks fractal dimension df using the relationship N ~ R d f . The 
fractal dimensions for the investigated threshold distributions r = —0.7 (o) and 
r = -0.5 (+) take the values df = 1.44 ± 0.10 and df = 1.39 ± 0.10 respectively. 
These values for df are consistent with the fractal dimension found in the two 
dimensional experiments on hydraulic fracturing of viscoelastic claysM 

The time record of the hydraulic pressure at the injection point is of con- 
siderable interest because on the one hand it is often experimentally measurable 
and on the other hand does it contain information about rheological and break- 
ing properties. If the considered crack is at rest the hydraulic pressure P will 
generally build up in time (loading). For an impermeable medium the pressure 
build up will correspond to the materials elastic response on the injected fluid 
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Figure 14: Temporal breaking sequence corresponding to the crack displayed 
in Fig. 11. The magnitude m(t) is defined as the number of per unit time step 



simultaneously broken beams. The simulation was stopped after 1500 time steps.. 
Note the clustering of breaking events and the large time intervals of quiescence.t3 



(and possible fluid flow effects within the crack which are neglect in the present 
case). We consider the fluid as incompressible implying a linear relationship be- 
tween crack volume and time, and due to the linearity and time independence 
of employed elastic equations a linear increase of pressure in time is the conse- 
quence. The associated slope is related to the materials mechanical stiffness (for 
the present crack). However, if the crack starts to move the stiffness decreases 
and so does the hydraulic pressure. There exists a cutoff for this pressure drop 
because the materials average cohesion is non-zero. Therefore the crack comes to 
rest again after some time period. Figure [H| displays such an 'erratic' pressure 
record. The strongest pressure fluctuations occur for small cracks (early times) 
at high frequency whereas lateron the loading periods become larger and larger. 
This is due to the rapidly decreasing stiffness as soon as the crack advances. The 
loading periods are periods of 'quiescence' (no breaking); the unloading periods 
correspond to 'bursts' (temporally localized groups of breaking events). This 
is exemplified in Fig. 14. The temporal correlations of such bursts have been 



analyzed in detail in Ref. |16 



5.2 The permeable medium 



This section considers hydro-fracturing of a permeable porous medium, where an 
incompressible fluid is injected at the center of a quadratic grid at a constant 
rate, qi n . The boundary conditions have been discussed in Sec. 2,4 . 

The dimensionless formulation of the pressure equation is first explained be- 
fore we show some examples. The considered pressure equation Eq. (^) is made 
dimensionless by considering the forces acting on the beams surrounding the 
injection cell at the center of the grid. Darcy's law yields the Darcy velocity 
v = —(k/fj)dp/dx in the x-direction. The Darcy velocity is also v in the y- 
direction when assuming symmetry. The forces acting on the beams surrounding 
the injection cell are F = Pddp/dx, because the pressure difference between two 
neighbor cells is dp = {dp/dx)l. The length of a beam is I and the width of a 
beam is d. The beam force F can be written F = pvl 2 d/k, when Darcy's law 
is used to replace the pressure gradient dp/dx with the Darcy velocity v. The 
volume flux into the center cell is initially qi n = vl 2 . The force F is written in 
terms of qi n rather than v , and we have F = pqi n d/k. This force is scaled with 
the unit force Fq = Ed 2 , (where E is Young's modulus), which can be written 
F/Fq = qi n /qo- The unit injection rate becomes qo = kEd/fi. 

This scaling shows two things. Firstly, it shows that a dimensionless injection 
rate of order one or larger implies strong enough forces on the center cell for 
deformations on grid to be noticeable. Secondly, the (dimensionless) average 
beam strength must be less than qin/qo for a crack to nucleate. 

The characteristic time becomes to = P-l 2 /{Ek), which will be as low as nano 
seconds for a grid with beam length I = m, permeability k ~ 10 -15 m 2 , 

Young's modulus E ~ 10 11 Pa, and fluid viscosity p, ~ 1CP 3 Pas. We consider a 
constant permeability throughout the material. 

The full pressure equation for an incompressible fluid Eq. (||) is not solved. 
The pressure equation is therefore not fully coupled to the volume changes in both 
the bulk and the fracture. However, it is still possible to couple the hydraulic 
pressure to the elastic deformations of the crack by considering conservation of 
fluid, as shown in Sec. 2.4. The (Laplacian) pressure equation is first solved for a 
unit pressure p* on the crack contour. The total Darcy flux D* is determined from 
this solution; the corresponding crack volume V* is computed from the associated 
elastic solution. The scaling parameter A (t) yielding the physical pressures and 



strains is computed each time step according Eq. (|il|), see Sec. 2A steps [fl|| 

The beams on the crack surface are checked after each time step to see if 
any of them are stretched beyond their strength. If there are any over-stretched 
beams, then the most over-stretched beam will be broken. We showed in the 
foregoing section that the advance of fracture is localized in time in form of 
bursts, separated by large time periods of quiescence (loading). The time-step 
used during the bursts is therefore "small" . On the other hand, the quiet periods 



Case 


Average beam strength 


Injection rate 


No 


(F th }/Fo 


qin/qo 


I 


1(T 2 


10" 3 


II 


1(T 3 


10~ 3 


III 


10" 4 


10~ 3 



Table 1: The numerical cases have the same injection rate qi n , but differ for the 
average beam strength (F t h)- 



between two bursts can be very long compared to the short time span of a burst. 
The time-step is therefore multiplied by a factor larger than 1 after time-steps 
where no beams were broken. This factor could be chosen so that ten such steps 
together yield a combined factor equal to for instance 2 or 10. Increasing each 
time-step between the bursts by a factor makes it possible to obtain a step length 
sufficiently long to get through these relatively long periods of pressure build-up. 
However, the time-step is immediately set to its minimum value once a new beam 
breaks. 

Three numerical examples are studied, which differ in the average beam 
strength (Fth)- The average beam strength in these cases and the injection rate 
qi n are given in the Table |]. All three cases have the grid size 100 x 100 and 
a beam strength distribution characterized by the exponent r = —0.9, compare 
Eq. 

A value for r close to —1 indicates that the beams have very broadly dis- 
tributed strengths. The employed material and geometry dependent constants 
are a = 1, b = 2.5 and c = 1200, when expressed dimensionless. (The scaling of 
these constants is such that a becomes 1.) 

The crack in case I, where qi n /qo "C (F t h)/Fo, is shown in Fig. 15 a. This 
figure shows the crack after it stopped growing. The hydraulic pressure is shown 
in Fig. 15b, where it is seen that the pressure reaches a stationary value. A small 
crack nucleated in this case, even though the average beam strength is 10 times 
larger than the injection rate. This is due to the strong heterogeneities (r = —0.9). 
The pressure follows an exponential rise in time according to the expression for 
the scaling parameter A in Eq. @. The first smooth rise of the pressure ends 
by breaking beams until no more beams are overstressed, and a second rise of 
pressure follows. This second rise leads to a stationary pressure. Notice that the 
first rise of pressure points towards a higher stationary pressure than the second 
rise. This is expected because the maximum possible pressure in a constant crack 
geometry, given a constant injection rate, is decreasing with increasing crack size. 
A large crack requires a lower stationary hydraulic pressure compared to a small 
crack in order to balance the injection rate to the hydraulic losses. The small crack 
seen in Fig. 15a is almost entirely due to the first burst. After the first burst, the 
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Figure 15: Case I has an average beam strength (Fth) /Fq = 10 -2 and an injection 
rate qi n /qo = 10 -3 . (a) (left) The final state of the branching crack structure and 
the corresponding Darcy flow field, (b) (right) The hydraulic pressure in the 
crack as a function of time 113 



pressure cannot reach a sufficiently high level any more to break more beams. The 
time r (see Eq. (|l0|)) is plotted in Fig. IS a. Because r measures the typical elastic 
response versus the typical flow response, it defines the characteristic time scale 
below which the (re)loading process appears to happen within an impermeable 
medium. It is seen that r is increasing in crack size, which is consistent with the 
pressure plot Fig. |l5|b. The characteristic time increases with the crack size even 
though the stationary pressure is decreasing. 

The case II shows crack propagation for an average beam strength equal to 
the injection rate, see Table [|. The forces from the fluid gradients on the beam 
structure are then comparable with the strength of the beams. The pressure 
during the cracking process is shown in Fig. |i~6|b, and it is seen that the pressure 
is rising almost linearly between the bursts of beam breaking. The small initial 
bursts were sufficient to generate a crack with a time constant r much larger than 
the time span shown in the pressure plot. This explains why the pressure rise 
appears almost linear rather than exponential in time as given by Eq. (^). The 
characteristic time r is shown in Fig. Eqb. The crack structure at the end of the 



time period in the pressure plot (Fig. [lqb) is shown in Fig. 16a. The crack in this 
case of broadly distributed beam strengths shows the growth of a branching crack 
structure. Note that within the crack 'fjords' the Darcy flow almost vanishes. 

Case III has an average beam strengths (F^/Fq one order of magnitude lower 
than the injection rate qi n /%- The injection rate is therefore high compared to 
strength of the structure. This is also seen from the pressure plot in Fig. 17b, 




Figure 16: Case II has an average beam strength (F^/Fq = 10~ 3 and an injection 
rate q% n /qo = 10~ 3 . (a) (left) The crack and the corresponding Darcv flow field, 
(b) (right) The hydraulic pressure in the crack as a function of time.0 



which is characterized by more frequent bursts compared to case II. Only a slight 
pressure increase is enough to trigger a new burst. The hydraulic pressure level 
of beam breaking in Fig. 17b is about one order of magnitude lower than the 
corresponding pressure level in figure Fig. 16 b. This is as expected from the 
average beam strengths of these cases. The characteristic time r for case III is 



shown in Fig. 18c, which also shows a large number of bursts. The crack in 



its final state is shown in Fig. |17|a. A branching crack structure is seen, where 
the branch close to the boundary seems to be growing faster than the other 
branches. This seems to be a boundary effect because the pressure gradients 
close to a boundary are larger. This effect favourizes growth of cracks close to 
the boundaries. 

The presented cases demonstrate several features of hydraulic fracturing in 
a permeable medium with broadly distributed cohesion, when a fluid is injected 
at a constant rate at the center of the test sample. They clearly show how the 
hydraulic pressure is increasing during "quiet" periods. The pressure increases 
according to Eq. The characteristic time r in Eq. ([R]) is computed, and 
it is seen that r becomes large compared to the time period between successive 
bursts. The pressure therefore appears to rise linearly between the bursts. The 
characteristic time r is also increasing with the crack size, which implies that the 
time periods between successive bursts get longer as the crack growths. This can 
also be seen from the slopes in Figs. |l5|b-(r^b which are decreasing with increas- 
ing crack size. The pressure drops during the bursts. Furthermore, the bursts 
also become more frequent when the average beam strength (F^/Fq gets low 




SCDD.C IDDDDO 



time [-] 



Figure 17: Case III has an average beam strength {Fth)/Fo = 1CP 4 and an 
injection rate qi n /qo = 10~ 3 . (a) (left) The crack and the corresponding flow 
field. Note the screening of the flow field within the crack branches, (b) (right) 
The hydraulic pressure in the crack as a function of time.O 



compared to the injection rate qin/qo- 



Conclusion 

We have reviewed some numerical results for hydraulic fracturing in porous ma- 
terials employing microstructural fracture growth models. The considered mod- 
els are not finite difference schemes of certain continuum mechanical equations. 
Rather they implement an approach to crack growth problems in the spirit of 
Monte Carlo methods of statistical mechanics. This allows to track quasistatic 
crack growth under very general conditions, i.e., for disordered systems on meso- 
scopic length scales. 
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Figure 18: The characteristic time r = J u* ■ d A/ / v* • <i A as a function of crack 
growth for the three cases of Table |l[ 
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